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Abstract 

We present new, efficient central schemes for multi-dimensional Hamilton- 
Jacobi equations. These non-oscillatory, non-staggered schemes are first- and 
second-order accurate and are designed to scale well with an increasing dimen- 
sion. Efficiency is obtained by carefully choosing the location of the evolution 
points and by using a one-dimensional projection step. First- and second-order 
accuracy is verified for a variety of multi-dimensional, convex and non-convex 
problems. 
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1 Introduction 

In this work we consider numerical approximations for solutions of multi-dimensional 
Hamilton-Jacobi (HJ) equations of the form 

+ H {x, <j>, V<ft) = 0, xeR N , (!•!) 

at 

subject to the initial data t = 0) — 

Hamilton-Jacobi equations are of special interest in a variety of applications, such 
as, e.p., optimal control theory, image processing, geometric optics, differential games 
and the calculus of variations. When the Hamiltonian does not depend on 0, solutions 
for (1.1) with smooth initial data will typically remain continuous but will develop 
discontinuous derivatives in finite time. Such solutions are not unique, and therefore 
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a mechanism is required for singling out a “physically relevant solution”, the viscosity 
solution. For convex Hamiltonians the viscosity solution coincides with the limit solution 
obtained by the vanishing viscosity method [11]. Extensions to general Hamiltonians 
were introduced by Crandall and Lions in [7], and have been systematically studied 
thereafter in a series of works [3, 5, 6. 25]. 

Hamilton-Jacobi equations are closely related to hyperbolic conservation laws. Yet 
while the literature on numerical methods for conservation laws is flourishing, very little 
attention is given to numerical methods for HJ equations. This is surprising given 
their increasing role in different applications. Crandall and Lions introduced in [8] first- 
order numerical approximations to the viscosity solution of a simplified version of (1.1), 
with a Hamiltonian that only depends on the derivative of 0. Discontinuous Galerkin 
(DG) methods for HJ equations were introduced in [10, 22]. Multi-dimensional DG 
schemes are based on transforming a scalar equation into a weakly hyperbolic system 
which is over- or under-determined, hence an additional least-squares step is required 
to single out a solution. High-order Godunov-type methods were introduced in Shu 
[32, 33], and were based on an Essentially Non- Oscillatory (ENO) reconstruction step 
that was evolved in time with a first-order monotone flux. The least dissipative flux, the 
Godunov flux, requires solving Riemann problems at cell interfaces. Central schemes 
avoid these difficulties by evolving the solution in smooth regions, i.e., by averaging 
over discontinuities. Such schemes have been widely studied for conservation laws, the 
prototype being the first-order Lax:- Friedrichs (LxF) scheme [9]. A one-dimensional 
second-order extension is due to Nessyahu and Tadmor [31]. Central schemes do not 
require Riemann solvers, which makes them suitable for solving systems of equations 
and for multi-dimensional problems. Extensions to two-space dimensions were done in 
[2, 14]; high-order central schemes were developed in [4, 23, 24, 28]; semi-discrete schemes 
that reduced the numerical dissipation and eliminated the staggering were developed in 
[16, 17, 19]. 

Central schemes have recently been extended to the HJ equations in [30], which 
applied the first- and second-order staggered central schemes of [14, 31] to HJ equations 
in one and two space dimensions. L 1 convergence of order one for this scheme was 
proved in [29]. In [18], a second-order semi-discrete scheme was presented, following the 
techniques for hyperbolic conservation laws [16, 19]. While less dissipative, this scheme 
requires the estimation of the local speed of propagation, which is computationally 
intensive in particular in multi-dimensional problems. In a later work, [17], the numerical 
viscosity was further reduced by computing more precise information about local speed 
of propagation. 

In this paper we derive non-staggered fully-discrete central schemes for approxi- 
mating solutions of (11). These methods combine the ideas of [18, 30], with several 
additional ingredients. Our scheme is presented as an W-dimensional algorithm which 
is designed with special consideration to performance and scaling to higher dimensions. 
We develop both first- and second-order accurate schemes. These schemes are based on 
a projection step which is one-dimensional regardless of the dimension of the problem. 
The methods described in this paper can also be thought of as the first step toward 
higher-order schemes, which is the subject of a forthcoming paper. 
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This paper is organized as follows: In Section 2 we develop our first and second 
order scheme in one dimension. Section 3 is the heart of the paper, where we gen- 
eralize these schemes to N dimensions, first introducing a multi-index notation, then 
deriving the location of the evolution points, and finally presenting the algorithm. Sec- 
tion 4 presents various examples, demonstrating the first- and second-order convergence 
of these schemes. We present a sample code implementation of our 2D second-order 
algorithm in the appendix. 

Acknowledgment: We would like to thank Ian Mitchell for helpful discussions and for 
suggesting the averaging strategy that is used to avoid staggering in the schemes. We 
also thank Volker Elling for helpful comments on early drafts of this paper. 


2 The One-Dimensional Scheme 

Consider the one dimensional Hamilton-Jacobi (HJ) equation 


<fi t + H(<f> x ) = 0, (2.1) 

subject to the initial data <p(x,t = 0) = (£>o(x). In order to approximate solutions of (2.1) 
we introduce a grid in space and time with mesh spacings, Ax and At, respectively. 
We denote the grid points by Xj = zAx and t n = nAt, and the fixed mesh ratio by 
A = At/Ax. Let denote the approximate value of 0(xj,t n ), and (</?*)" denote the 
approximate value of the derivative 4> x (x,, f n ). We define (A<p)" + i := 

Given <p™, an approximate solution at time f n , the approximate solution at the next 
time step f n+1 , is obtained as follows: 

1. Reconstruct a continuous piecewise-polynomial from the data, and sample it 
at the half-integer points, {xi+ 1 / 2 }, to obtain the values of <p” +i and its derivative, 

(<Ar)"- The order of the polynomial is related to the overall order of accuracy of 
the method. 

2. Evolve <p n 1 by solving (2.1) from time t n to time £ n+1 , obtaining ■ This 

t+ 2 ^ y~2 

evolution is done at the half-integer grid points where the reconstruction is smooth, 
so long as the CFL condition A | H' (y? x )| < 1/2 is satisfied. 

3. Project y? n+ i 1 back onto the integer grid points {x*} to get y>” +1 . 


2.1 A First Order Method 


The derivation of the first-order method starts by reconstructing a piecewise-linear in- 
terpolant of the form 



Ax 


(x - X.) 




<p(x,t n ) := 


(2.2) 
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where Xi+iM is the characteristic function of the interval [xj,x i+ i)- The values of the 
interpolant (2.2) and its derivative at the half-integer grid points, x i± i are 


<p?±i = ^ ± o ( A t?)r± 


(A^i 

MU = 


Ax 


Integrating (2.1) in time from t n to t, n+1 at x i± ^ gives 

<4 = - AtH (mu) = ± \ ( A »>>"±) - A tH ( fe ^ i ) - 

Finally, we project the evolved solution back onto the original grid points. For a first- 
order method it is sufficient to average V&i/v 



The intermediate values Vi+ 1/2 are t ^ ie same 35 those computed in the first-order 
method in [31], so in one dimension we only add the projection step. This eliminates the 
grid staggering in [31] with little computational cost since no additional flux evaluations 
are required. 


2.2 A Second Order Method 

The second-order scheme is based on a piecewise-quadratic interpolant of the form 


V 


(x,t n ) := 


Vi + 


(A^)”, 

Ax 


(x — Xi ) + 


lp(Art" + j 

2 (Ax) 2 


(x - x,) (x - x, + i) 


Xi+l-(2.4) 


Here, T> is a limiter whose goal is to prevent oscillations while maintaining the order of 
accuracy of the method. There are various possibilities for choosing such limiters (see 
[37]). One such example is the Min-Mod limiter, 


Vfi := MM 


»(/.+! - Si), \(f ,+ 1 - «(/.+! - Si) 


1 < e < 2, 


where the Min-Mod function is defined as 

{ minj {x_j} , if all Xj > 0. 
maxj{xj}, if all Xj < 0, 

0, otherwise. 

Since the solution of the HJ equations generally has a discontinuous first-derivative, we 
follow [18] by limiting the second-derivative. With the Min-Mod limiter, the second 
derivative is approximated by D(A<^)" +1 , 2 /(Ax) 2 , where 

K(Art” + i = MM [fl ((Art" +§ - (Art?*}) , \ ((Art” +f - (Art".,) , 

6 ((Art" + , - (Art"-i)] ■ 
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Sampling (2.4) and its derivative at the half-integer grid points gives 


1 1 (A^)Li 

± 5 (A*)r ± i - gO (A^)” ± 1 . (9.)U = 


Ax 


We integrate (2.1) from time t n to time t n+l using the second-order midpoint quadra- 


ture 


J H (*>, (x, ± i ,(“)) it ~~ A tH (fe)"4) • 

The required mid- values, tp x ^x i± i , t n+ 5 J , can be predicted using a Taylor expansion, 
‘Px (^i±i^ n+ ^) =<Px + ^Atifit x (x i± !_ + O ((At) 2 ) = 

= <p x (<Px (^±1,*")) ¥>XX (^i±I,< n ) + O ((At) 2 ) S3 


(A<P)r ± i i /(a^ ± a p(A^r ± r 

Ax 2 AjW l Ax / Ax 


which leads to 


*2* = <&* - 


/ (Ay?)" 

Ax 


‘±5 


lwr/ /(ACi^M ± |' 

2 XH \ Ax j Ax 


(2.5) 


Finally, we project (2.5) back onto the integer grid points using a quadratic interpolant 




where (Ay?)" +1 = <p n +} — <p n +i , and 
2+ 2 2 

© (A v )r 1 = mm [e ((A v )"« - (Av)r +t ) , 1 ((/Mm - (A?)”* 1 ) , 


» «a^ +i - (Aic, 1 ) \ 


Remark. We would like to note that even in the ID scheme, there are several differences 
between our method and the second-order scheme in [30]. A second-order interpolant is 
used to re-project the evolved fields back onto the original grid points, resulting with a 
non-staggered grid compared with the staggered scheme in [30]. Also, we follow [18] by 
applying the nonlinear slope limiters to the second derivative. 
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3 Generalization to N Dimensions 

We are concerned with approximating solutions of the iV-dimensional HJ equation of 
the form 

4>t + H{V<t>) = 0, xeR N , (3.1) 

subject to the initial data <t>(x, t) = 4>o(x). 

In Section 3.1 we introduce a multi-index notation, which allows a presentation 
that nicely parallels the one-dimensional case. We then compute the optimal location 
of the evolution points. Sections 3.2 and 3.3 develop the first- and second-order N- 
dimensional schemes. The first-order method in §3.2 below applies as is to the case 
where the Hamiltonian H depends also on x and 0. We extend the second-order method 
of §3.3 to this more general case in a remark. 

3.1 Preliminaries 

A multi-index notation. We define the multi-index a; = (au, ct 2 , • • • , (*n), and denote 
by x a , the point x a — ^Xq,\ XaJ, • • ■ , £ M^. Here x^ denotes the fc-th coordinate 

of x so XaJ = a^Ax^K For example, in the conventional three-dimensional notation 
with indices i, j and k and components (x, y , z), a = ( i,j , k ) and x a = (x*, z*. ). 

For a given a we define the special multi-indices a ± e* := (c*i, . . . , Qjt ± 1, . . . , ), 

which denotes an increment in the k direction. Then ip'l = y>{xa} ■, x^l , - . - ,x^J; t n ) and 
t Pa±e k = ^(xi 1 /, . . . , x l £ ± Ax( fc \ . . . , x^J; t n ). Finally, we denote the evolution points 
with the multi-indices ± := (cq ± a, Q 2 ± 0 , . . . , ajv ± a), for some constant a , so that 
‘Pdt — <p{%ai i aAx^\ Xal ± aAx^ 2 \ . . . , x^J ± aAx^\ t n ). 

The location of the evolution points. We would like to determine the optimal 
location of the evolution points. For simplicity we assume a grid point at the origin, 
x = (0, 0, . . . , 0), and scale the coordinates such that VA;, Ax^ = 1. The two evolution 
points will then be located at x± = (±a, ±a, . . . , ± 0 ) for a constant a that is yet to be 
determined. 

Consider the evolution point x+, which is at a distance y/Na from the origin. The 
value of ip at this point will be based on a polynomial that is constructed inside the 
(hyper-) volume bounded by the coordinate planes and the (hyper-) plane x ^ = 1- 

There will be discontinuities in the first-derivative of the piecewise-polynomial inter- 
polant tp(x, t n ) along the sides of this hyper- volume. Since we evolve the solution in 
smooth regions, an optimal choice of the evolution points is at an equidistant location 
from these boundaries. The diagonal line (s, s, . . . , s) (for some parameter s) intersects 
the (hyper-)plane YliLi x ^ = 1 at s = l/N, or at the point x p = (-^, which 

is at a distance l/y/N from the origin (see Figure 3.1). 

The optimal choice is to require that the evolution points are equidistant from the 
coordinate planes and the intersection point x p . The distance from x + to all the coordi- 



Central Schemes for HJ Equations 


7 


nate planes is a. The distance from x+ to x p is 1/y/N-ayfN, therefore the requirement 
that x + be equidistant from the coordinate planes and x p is 

1 

a = t=. 

N + s/N 

The evolution points in [30] were chosen as a = 1/4, which places them equidistant 
between the origin and the intersection point x p . In N dimensions this choice generalizes 
to a = (2 N)~ 1 . In our case, when N — 2, a = (2 + V^) -1 ~ 0.29 which is about 15% 
larger than the choice a = 1/4. When N — 3, a = (3 + V^) -1 « 0.21 which is about 
30% larger than the choice a = 1/6. Thus the optimal choice of a will allow larger mesh 
ratios, leading to larger time steps and less dissipation. 



Figure 3.1: The location of the evolution points x± and x p in 2D. 


3.2 A First Order Method 

For simplicity we assume that the spacing is identical in every direction, i.e., Ax (fc) = Ax, 
for all k. Generalization of the methods below to the case where ^ Ax^ for 

k ^ j is straightforward. We define the forward- and backward-differences in the k - th 
component as A := -y?" and A := ~Va-e k i respectively. At each grid 

point x a we reconstruct two linear interpolants that are valid in the two hyper-quadrants 
that contain the points x± = x a ± (a, . . . , a) Ax, 

v± (*. n ■■=<&+'£, ^ ■ <X2) 

k = l 

In order to compute the solution at the next time step at x a , we first compute the 
solution at time t n+l at the evolution points x±, and then average these two values. The 
value of the linear interpolant (3.2) at x± is 

N 

rix±,n = l fii±aY,*t'Pi. 

1 
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and its derivative is 


(V< := 


Afy; 

Ax ’ Ax 


Hence, the values at the evolution points x± at the next time step, t n ^ x , are given by 

t n + l 

(p{x±,t n+1 ) = <p(x±,t n ) - / H (Vip(x±,t n ))dt « p(x±,t n ) - AtH ((Vcp)i) = 

Jt n 


= *»±a 


Ajv; a^; 


The value at t n+l at x a is finally obtained by averaging y?± +1 := <^(x±,t n+1 ) (compare 
with (2.3)), 


¥>" +1 = 2^" +1 + ^- +1 ) = 


\/t=l fc=l / 


At K< \ , 

2 Ai Ax i V Ax ’ ' ’ Ax 


3.3 A Second Order Method 

For simplicity we assume again that the mesh spacing is identical in every spatial di- 
rection, i.e. , Ax^ = Ax, for all k. Similarly to the ID case in Section 2.2, the N- 
dimensional second-order method is based on a piecewise-quadratic polynomial. For 
every grid node we reconstruct two iV-dimensional quadratic interpolants: '~p + (x, t n ) 
for the hyper-quadrant along the positive diagonal, and <p_ (x, t n ) along the negative 
diagonal (see Figure 3.1), 


*>*(*.« ") :=rf + E x«‘>) + 


+\ £ 70 s (*“> - 4‘>) (* w - 41.) + 

1 N N A±,,n 

+ l - y y (x (j) - xg>) (* (fc) - ®«) . 

2 ^tT (Ax ) 2 v } 

The Min-Mod limiter in the j-th direction acting on A ± <£>” is 

t>a i< = mm [e (a - A±<) , 

5 (a?*S«, - Af^) ,«(Af^ - A±^_,.)] . 


(x<‘> - x?>) (x<‘> - x<2..) + 
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so that (Ax) 2 approximates the second derivative d 2 <p (x a ,t n ) j dx^dx^ . 

Now x± } - ia ] = iaAi, x+* - xi+ et = (a - 1) Ax and xf fe) - x^l ek = - (a - 1) Ax 
so evaluating (3.4) at x± 


< P ± '■= < P ±{ x ±, t n ) 


= & ± a £ zta + y, j\a + f E £ 


fc=l 


fc=l 


7 = 1 ^=1 


The approximation to the first derivative of (3.4) in the p-th direction is given by 


dipt(x,n _ AjvS PpA?»5 [/>) ,„A I' M (p) 'll 

“ at + x ° ) + \ z± x ^)\ 


+ 


+ 


2 (Ax) 


(-2 £ [D p Afrf + X>*AjvS] (xL‘* - *?>) . 




which when evaluated at x± is 
f d<p = dy± (j, C) 




± 

AJ ± 2a-2 PpA^; ± a A DpAf 

^ fcssl 
*7*p 


<9xtp) 


*± 

± 

p 


Ax 


Ax 


Ax 


The approximation to the second derivative is given by 

( ay± y = p,a*< + p,a^; 

\dx^dx^ J 2 (Ax) 2 

The solution at the next time step at the evolution points <p± +1 i s obtained by 
evolving the reconstruction (3.4) according to (3.1). The integral of the Hamiltonian 

is approximated by a second-order midpoint quadrature, /; H(v<p(x i± ,n)dt « 


Atff 


which at the evolution points gives 


(<Vvp)l +i ), 
v4 +1 =<-Atff ((V*>)i +i ). 


(3.5) 


Here ( V 7 ?)" := (( 5707 )^. ■ ■ ■ . ( 77 ^ 7 ) ±j denotes the approximation to the gradient at 
x±. The mid-values in time can be estimated via the Taylor expansion using (3.1), 


dp 

dx^ 




dp 

dx^ 


dX) (X± ’ + 2 dxMdt (X± ’ <n) + ° ^ 

N 


(3.6) 


<**■ r » gjtim (*±-n+o (a^) . 


2 
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Hence, 



Ax 


2a - 1 V p A^(pl 
2 Ax 


, a^V p A±y£ + £> fc Aj^ 
aJ 

fc=l 

fc^p 


At 

T 


E 


°dx£ 


VpA^ + VkAp^_ 
2 (Ax) 2 



Figure 3.2: The location of the points x+, x_, X(_i) + and X(+i)_ used in the projection 
step along with the distances Do, D + and D _ in the two dimensional case. 

All that remains is to project (3.5) back onto the original grid points, x Q . This 
projection is one-dimensional regardless of N. We use the four evolution points x + , x_, 
x ( _i)+ := (xi^ - Ax (1) + aAx (1) , xi 2) - Ax (2) + aAx (2) , . . . ,xi N) - Ax (N) + aAx (lV) j and 

x (+1) _ := (xL 1} + Ax (1) - aAx (1) , xL 2) + Ax (2) - aAx (2) , . . . , x ( a N) + Ax (JV) - aAx^ (see 

Figure 3.2). The distances between the evolution points are D 0 |x+-x_| = 
2a\fN Ax, D + := |x (+1 )_-x + | = (1 - 2a) \fNAx, and D- := |x(_i)+-x_| = D+. 
We then define the approximations to the first derivative along the diagonal, 

wr 1 , = <■' - ^ w: +l <i‘- - 

Do Do D+ D+ 

(A 

(“¥>)_ _ ¥4 ^(-1)+ 

D_ £>_ 
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The approximation to the second derivative is the limited difference V(dp) 0 / D 2 , 

1 

2 


where D 2 = | (Do + D + ), and 


V(dp)^ + = MM 

6 


9 


(d<p) 


n+1 


(<v>r‘ \ 1 1 w 


n+1 


\n+l ' 


D + Do }'2\ D + 

(dip) r 1 w +l ' 

D 0 D_ 


D_ 


(3.7) 


The approximated value at the next time step t n+1 at the grid point x a is therefore 
given by 


<p n a +X = pl +1 + — 7^ (x 0 -x_) + :P ^ 0 (x a -x_)(x a -x+) 


\n+l 


9_ +1 + 


D 0 2 D 2 

yf ~ V >- +1 A) D(d^)” +1 ^o 


Do 


2D, 


where (p± +1 is given by (3.5) and V (dp) q +1 is given by (3.8). 

Remarks. 

1. If the Hamiltonian H depends also on x and <£, then (3.5) becomes 
Pl +l =p n ± - AtH (x,pTK (V<^) , 


where 


At 


<P± 5 = <P± ~ ~j H ( f ’ ^±> ( V< ^)±) > 


and the Taylor expansion (3.6) contains the additional term 


At 

~2 




(i±,t n ) 


2. We would like to stress that the fully-discrete scheme in [18] that was derived as an 
intermediate step in developing the semi-discrete scheme, was only first-order in 
time. Moreover, in 2D our scheme is based only on two flux evaluations compared 
with four flux evaluations in [17, 18]. It also does not require any estimation of 
the local speed of propagation at every grid point (as required in [17, 18]) at the 
price of being more dissipative. 

We would like to summarize the second-order iV-dimensional algorithm for a general 
Hamiltonian H(x, 4>, V0) : 
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Algorithm 3.1 Let the distance of the evolution points from the origin be a — 

1. For each grid node x a and each k compute 
KfI = Vl+e k ~ Fa and A*</£ = </£ - 

2. For each grid node x a and for each j and k compute 


T>Ai< = mm [» (aJ^ +(j 


-xta 


~ i&tFl+ej 


e (a to - Atvl-.,)] ■ 


3. For each grid node x Q compute 




a (a — 1) 




j = 1 


and for each p compute 

( d(p \ n _ Ap 2a — 1 £pAp <Pg ± a Dp A^p” + DkApp^ 

\dxbl) ± ~ Ax 2 Ax 2 Ax 

k*p 


(_ dp _ y +h 

\c>x^ / ± 


(£k)\~T + 

+ ±H( f,^,(V v )i) (J^)” + 
z^ d ^_ y ^ [ 2 (Ax) 2 


F n a t l = Fl- AtH ((Vp)l H ) , 

where H ((V</?)±) = H ((^y)” , • • • , (a^o)2) • 

4. Let D 0 = 2ay/N Ax, D+ = D_ = (1 - 2a) \/N Ax. For each x Q compute 

(<« +1 wS 1 - mV 1 _ <+«- - v"X' _ <+ ~ <-V 

D 0 Do ’ D+ D + £>_ D _ 
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(where a ± 1 is the multi-index (aq ± 1, . . . , aw ± 1) ) 


®(<V)o" + ‘ = 


MM 


9 


w ; +1 

D + 


m q +1 > \ 

Do ) 


9 


W o +1 

Do 


W2 +1 \ 

) 


l/«‘ (<V): +1 \ 

2 D + D. )' 


5. For each x a compute 


<P, 


” +1 - ^ KV + <*-) 




4 Numerical Examples 

We demonstrate the schemes developed in Section 2 and Section 3 with several examples. 
Most of these examples are standard test cases that can be found, e.g., in [18, 30, 33]. 


Example 1: A convex Hamiltonian 

We start by testing the performance of our schemes on a convex Hamiltonian. We 
approximate solutions of the one-dimensional equation 

(fit + - ((fi x + l) 2 = 0, (4.1) 

subject to the initial data (fi(x, 0) = — cos(7rx) and to periodic boundary conditions on 
[0, 2]. The change of variables, u (x, t) = (fi x (x, t) + 1, transforms the equation into the 
Burger’s equation, u t - 1-| ( u 2 ) x = 0, which can be solved via the method of characteristics 
[33]. As is well known, Burger’s equation generally develops discontinuous solutions even 
with smooth initial data, and hence we expect the solutions of (4.1) to have discontinuous 
derivatives. In our case, the solution develops a singularity at time t — n~ 2 . 

The results of our simulations are shown in Figure 4.1 and Figure 4.2. The order 
of accuracy of these methods is determined from the Z^-norm of the error (see [29]). 
The results before the singularity, at T = 0.8/7T 2 , are given in Table 4.1, and after the 
singularity at T = 1.5/ 7r 2 in Table 4.2. 

In 2D we solve a similar problem 

(fit + 2 (0x + (fi y + l) 2 = 0, (4-2) 

which can be reduced to a one-dimensional problem via the coordinate transformation 

The results of the second-order calculations for the 

initial data cfi(x,y, 0) = — cos (7r(x + y)/2) = — cos(7rf) are shown in Figures 4. 3-4. 4. 
The convergence rates for the first- and second-order 2D schemes before and after the 
development of the singularity are shown in Tables 4. 3-4. 4. 
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N 

First-order method 

Second-order method 

L l -error 

L l -order 

L l -error 

L 1 -order 

40 

0.0342 

- 

6.14xl0 -3 

- 

80 

0.0172 

1.00 

1.53X10" 3 

2.00 

160 

0.0082 

1.06 

3.62xl0~ 4 

2.08 

320 

0.0040 

1.03 

8.77xl0- 5 1 

2.04 


Table 4.1: ZJ-errors for the ID convex HJ problem (4.1) before the singularity formation. 
T = 0.8/tt 2 . 


N 

First- order method 

Second-order method 

L l -error 

L l -order 

L l -error 

L l -order 

40 

0.0788 

- 

0.0110 

- 

80 

0.0379 

1.06 

2.70x 10 -3 

2.03 

160 

0.0183 

1.05 

6.79xl0" 4 

1.99 

320 

0.0091 

1.00 

2.28x10-^ 

1.57 


Table 4.2: ZJ-errors for the ID convex H J problem (4.1) after the formation of the singularity. 

T= 1.5/tt 2 . 


N 

First-order method 

Second- order method 

L 1 -error 

L l -order 

L l -error 

L l -order 

50 

0.299819 

- 

0.0261816 

- 

100 

0.144967 

1.05 

0.0061162 

2.10 

200 

0.0712024 

1.03 

0.00146501 

2.06 


Table 4.3: ZZ-errors for the 2D convex HJ problem (4.2) before the singularity formation. 
T = 0.8/tt 2 . 


N 

First-order method 

Second-order method 

L l -error 

L l -order 

L l -error 

L l -order 

50 

0.44437 

- 

0.0356362 

- 

100 

0.215055 

1.05 

0.0101615 

1.81 

200 

0.104421 

1.02 

0.00237225 

2.10 


Table 4.4: /^-errors for the 2D convex HJ problem (4.2) after the singularity formation. 
T = 1.5/tt 2 . 
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We proceed with a 3D generalization of (4.2), 

0t + ^ (^z + <Py +02 + 1) = 0) (4-3) 

subject to the initial data <p ( x , y, 0) = — cos (rc(x + y + z)f 3). The convergence results 
for the first- and second-order 3D schemes before and after the singularity formation 
are given in Tables 4. 5-4. 6. 


N 

First-order method 

Second-order method 

L 1 -error 

L l -order 

L l -error 

L l -order 

50 

5.932 

- 

0.672 

- 

100 

2.838 

1.06 

0.155 

2.12 j 

200 

1.76 

0.69 

0.041 

1.92 


Table 4.5: L l -e rrors for the 3D convex HJ problem (4.3) before the singularity formation. 
T = 0.08. 


N 

First-order method 

Second-order method 

L l -error 

L l -order 

L l -error 

L l -order 

50 

8.801 

- 

0.776 

- 

100 

4.148 

1.09 

0.171 

2.18 

200 

2.138 

0.96 

0.055 

1.65 


Table 4.6: ZJ-errors for the 3D convex HJ problem (4.3) after the singularity formation. 
T = 0.152. 


Example 2: A non-convex Hamiltonian 

In this example we deal with non-convex Hamilton- Jacobi equations. In ID we solve 

(pt - cos (<p x + 1) = 0, (4.4) 

subject to the initial data, <p(x, 0) = — cos(7nr), and periodic boundary conditions 
on [0,2]. In this case, (4.4), has a smooth solution for t < 1.049/7T 2 , after which a 
singularity forms. A second singularity forms at f ~ 1.29/7T 2 . The results are shown in 
Figures 4. 5-4. 6. The convergence results before and after the singularity formation axe 
given in Tables 4. 7-4. 8. 

In 2D we solve 

(pt - cos (<p x + <p y + 1) = 0, (4.5) 

subject to the initial data, (p ( x , y, 0) = — cos (7r(x + y)/ 2), and periodic boundary con- 
ditions. The results are shown in Figures 4. 7-4. 8. The convergence results for the first- 
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N 

First- order method 

Second-order method 

L l -error 

L l -order 

L l -error 

L l -order 

40 

0.0406 

- 

7.71 x 10 _a 

- 

80 

0.0206 

0.99 

2.00 xbH 

1.95 

160 

0.0099 

1.06 

4.75 xKT 4 

2.07 

320 

0.0048 

1.03 

1.16xl0" 4 

2.04 


Table 4.7: L 1 -errors for the ID non-convex HJ problem (4.4) before the singularity formation. 

T = 0.8/tt 2 . 


N 

First-order method 

Second-order method 

L l -error 

L l -order 

L l -error 

L l -order 

40 

0.1057 

- 

0.0248 

- 

80 

0.0515 

1.04 

6.92 x 10 -2 

1.84 

160 

0.0248 

1.05 

1.89xl0 _a 

1.87 

320 

0.0121 

1.04 

4.96 xl0~ 4 

1.93 


Table 4.8: Z^-errors for the ID non-convex HJ problem (4.4) after the singularity formation. 
T = 1.5/7 r 2 . 


N 

First- order method 

Second-order method 

L l -error 

L l -order 

L l -error 

L l -order 

50 

0.4046612 

- 

0.0402387 

- 

100 

0.195242 

1.05 

0.00999105 

2.01 

200 

0.0958819 

1.03 

0.0024644 

2.02 


Table 4.9: L^errors for the 2D non-convex HJ problem (4.5) before the singularity formation. 
T = 0.8/tt 2 . 


N 

First-order method 

Second-order method 

L 1 -error 

L l -order 

L l -error 

L l -order 

50 

0.705644 

- 

0.0837231 

- 

100 

0.334981 

1.07 

0.0215121 

1.96 

200 

0.161502 

1.05 

0.00555327 

1.95 


Table 4.10: L l -e rrors for the 2D non-convex HJ problem (4.5) after the singularity formation. 
T = 1.5/tt 2 . 
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and second-order 2D schemes before and after the singularity formation are given in 
Tables 4.9-4.10 and confirm the expected order of accuracy of our methods. 

The extension of (4.5) to 3D reads 

0 t - cos {<p x + <t> y + (f> z + 1) = 0, (4.6) 

The initial data is taken as <f> (x, y, 0) = — cos (ir(x + y + z)/ 3) . The convergence rates 
for the first- and second-order 3D schemes are given in Tables 4.11-4.12. 


N 

First-order method 

Second-order method 

1 '} -error 

L l -order 

L l -error 

L l -order 

50 

5.071 

- 

0.603 

- 

100 

2.55 

0.99 

0.150 

2.01 

200 

1.142 

1.16 

0.035 

2.10 


Table 4.11: Zd-errors for the 3D non-convex HJ problem (4.6) before the singularity forma- 
tion. T — 0.08. 


N 

First-order method 

Second-order method 

L l -error 

L l -order 

L l -error 

L l -order 

50 

8.594 

- 

1.100 

- 

100 

4.26 

1.01 

0.299 

1.88 

200 

1.913 

1.16 

0.075 

1.99 


Table 4.12: Zd-errors for the 3D non-convex HJ problem (4.6) after the singularity formation. 
T = 0.152. 


Example 3: A linear advection equation 


In this example we solve the ID linear advection equation, i.e., the Hamiltonian is taken 
as H (4> x ) = <f> x . We assume periodic boundary conditions on [—1, 1], and take the initial 
data as 4> (x, 0) = g (x — 0.5) on [—1,1], where 


<?(*) = 


V3 9 2t 

T + 2 + T 


(x + 1) + h(x), 


h(x) 


/ 


l 


2 cos (y x 2 ) — \/3, 

| + 3 cos (2nx ) , 
y — 3 cos ( 2nx ) , 

| (28 + 47T + cos (37rx)) - cos (fx) , 


— 1 < X < — 

— | < x < 0, 

0 < x < 

1 <x< 1. 


This example was designed as to be similar to the one used in [12]. An error in that 
reference prevents us from repeating it exactly. The results of the second-order method 
are shown in Figure 4.9. The dissipation effects axe visible in the round-off of the corners. 
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Example 4: 2D eikonal equation in geometric optics 

We demonstrate the results obtained with the 2D scheme on the non-convex problem 


( 4>t+ V^I + + 1 — 0’ (4.7) 

\ 0 (x, y, 0) = \ (cos (2ttx) - 1) (cos (2 ny) - 1) - 1. 

This model arises in geometric optics [15]. The results of our second-order method at 
time T = 0.6 are shown in Figure 4.10, where we see the sharp corners that develop in 
this problem, in agreement with the results in [30]. 


Example 5: Optimal control 

We solve a 2D problem with a more general Hamiltonian of the form H(x,y,\7<fi). This 
is an optimal control problem related to cost determination [33]. 

{ 4> t - sin (y) 4> x + sin (x) 0 y + \4> y \ - \ sin 2 (y) - 1 + cos (x) = 0, / 4 g\ 

\ 0(x,y, 0) = 0. 

This example develops a complex singularity structure. The result of our second- 
order scheme is in qualitative agreement with [30] as can be seen in Figure 4.11. 



Figure 4.1: Example 1. The ID convex Hamiltonian (4.1) before the formation of singu- 
larities. T = 0.8/ 7 r 2 . N = 100. Shown are the first-order approximation, the second-order 
approximation and the exact solution. 
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0 
0 
0 
0 

' 0 . 

- 0 . 

‘ 0 . 

- 0 . 

■ 1 . 

Figure 4.2: Example 1 . The ID convex Hamiltonian (4.1) after the formation of singular- 
ities. T = 1.5/7T 2 . N = 100. Shown are the first-order approximation, the second-order 
approximation and the exact solution. 




Figure 4.3: Example 1. The 2D convex Hamiltonian (4.2) before the formation of singular- 
ities. T = 0.8/ 7T 2 . N — 40 x 40. 
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Figure 4.4: Example 1. The 2 D convex Hamiltonian (4.2) after the formation of singularities. 
T = 1.5/tt 2 . N = 40 x 40. 



Figure 4.5: Example 2. The ID non-convex Hamiltonian (4.4) before the formation of 
singularities. T = O. 8 / 7 r 2 . N — 100. Shown are the first-order approximation, the second- 
order approximation and the exact solution. 
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Figure 4.6: Example 2. The ID non-convex Hamiltonian (4.4) after the formation of 
singularities. T = 1.5/tt 2 . N = 100. Shown are the first-order approximation, the second- 
order approximation and the exact solution. 



Figure 4.7: Example 2. The 2D non-convex Hamiltonian (4.2) before the formation of 
singularities. T = 0.8/ir 2 . N = 40 x 40. 
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Figure 4.8: Example 2. The 2D non-convex Hamiltonian (4.2) after the formation of 
singularities. T = 1.5/7 r 2 . N = 40 x 40. 




Figure 4.9: Example 3. A ID linear advection problem. N = 400. 
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Figure 4.10: Example 4. The 2D eikonal equation (4.7). N = 40 x 40. Left: the initial 
data. Right: the solution at T = 0.6. 



Figure 4.11: Example 5. The 2D optimal control problem (4.8). T = 1. N = 40 x 40. 
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Appendix A: A Second-Order 2D MATLAB Imple- 
mentation 

We present an example of a MATLAB implementation of Algorithm 3.1 in 2D. The 
following code provides the core computations for a single time step, but does not include 
any handling of boundary conditions. The functions H(phix, phiy), HphixCphix, 
phiy) and Hphiy(phix, phiy) return, respectively, the Hamiltonian H (tp x ,<p y ) and 
the derivatives of the Hamiltonian and ■ The function minmod is the 

slope limiter defined in Section 2.2. This code is somewhat optimized for speed at the 
expense of storage, and h := Ax, k At. 


al « 1/(2 + sqrt (2) ) ; 
a2 - al*(al-l)/2; 
a3 ■ al*al/2; 
a4 - (2*al-l)/2; 
a5 - al/2; 
h2 - h*h; 
twoh2 - 2*h2; 
kover2 » k/2; 
for i*l:N 

for j-l:N compute the first differences 
dphix(i,j) * phi(i+l,j) - phi(i,j); 
dphiy(i.j) * phi(i,j+l) - phi(i,j); 

end 

end 

for i*l:N 

im * i-1; 

ip * i+1; 

for j»l:N */, compute the limited second differences 
jm - j-1; 
jp - j+i; 

dphixx(i,j) - minmod (dphix( ip, j) - dphix(i,j), dphix(ip.j) - dphix(im,j), dphix(i,j) 

dphiyx(i,j) ■ minmod (dphiy( ip, j) - dphiy (i,j), dphiy ( ip, j) - dphiy(im.j), dphiy(i,j) 

dphixy(i , j ) * minmod(dphix(i , jp) - dphix(i,j), dphix(i,jp) - dphix(i,jm), dphix(i.j) 

dphiyy(i,j) * minmod (dphiy (i ,jp) - dphiy(i.j), dphiy(i,jp) - dphiy(i,jm), dphiy(i,j) 

end 


dphixCim, j)) ; 
dphiy (im, j)) ; 
dphix(i , jm) ) ; 
dphiy(i, jm)) ; 


end 

for i-l:N 

for j*l:N 

im - i-1; 

jm - j-1; 

dphixyp - dphiyx(i,j) + dphixy ( i , j ) ; 
dphixym ■ dphiyx(i,jm) + dphixy(im, j) ; 

phip - phi(i,J) + al*(dphix(i, j) + dphiy(i.j)) + a2*(dphixx(i, j) + dphiyy(i,j)) + a3*dphixyp; 
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phim - phi(i,j) - al* (dphix(im, j) + dphiy(i.jm)) + a2*(dphixx(im, j) + dphiyy(i , jm) ) + a3*dphixym; 

phixp - (dphix(i.j) + a4*dphixx(i , j) + a5*dphixyp) /h ; 

phiyp * (dphiy(i,j) + a4*dphiyy(i , j) + a5*dphixyp) /h; 

phixm ■ (dphix(im.j) - a4*dphixx(im , j) - a5*dphixym)/h; 

phiym - (dphiy(i,jm) - a4*dphiyy(i , jm) - a5*dphixym)/h; 

phixtp ■ -Hphix (phixp, phiyp) *dphixx(i, j)/h2 - Hphiy (phixp, phiyp) *dphixyp/twoh2 ; 
phiytp ■ -Hphix (phixp, phiyp) *dphixyp/tvoh2 - Hphiy (phixp, phiyp) *dphiyy(i, j)/h2; 
phixtm - -Hphix (phixm, phiym) *dphixx(im,j)/h2 - Hphiy (phixm, phiym) *dphixym/twoh2; 
phiytm - -Hphix (phixm, phiym) *dphixym/tvoh2 - Hphiy (phixm, phiym) *dphiyy(i ,jm)/h2 ; 
phihatp(i.j) » phip ~ k*H (phixp ♦ kover2*phixtp, phiyp + kover2*phiytp) ; 
phihatm(i,j) - phim - k*H (phixm + kover2*phixtm, phiym + kover2*phiytm) ; 

end 


end 

Dp * (l-2*al)*h^sqrt (2) ; 

Dm - Dp; 

D2 - (DO + Dp)/2; 

D02o8D2 * D0*D0/(8*D2); 
for i*l:M 

for j-l:N 

ip ■ i+1; 

jp - j+i; 

im - i-1 ; 
jm - j-t; 

dphiO ■ (phihatp(i , j) - phihatm(i , j) )/D0; 

dphip ■ (phihatm(ip, jp) - phihatp(i , j) )/Dp; 

dphim - (phihatmCi , j ) - phihatp(im, jm))/Dm; 

d2phi ■ minmod (dphip - dphiO, dphip - dphim, dphiO - dphim); 

phi(i,j) - (phihatp(i, j) ♦ phihatm(i , j ) )/2 - D02o8D2*d2phi ; 

end 

end 
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